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Abstract. We present limits for the compactification scale in the theory of Large Extra Dimensions 
(LED) proposed by Arkani-Hamcd, Dimopoulos, and Dvali. We use 11 months of data from the 
Fermi Large Area Telescope (Fermi-LAT) to set gamma ray flux limits for 6 gamma-ray faint neutron 
stars (NS). To set limits on LED we use the model of Hannestad and Raffelt (HR) that calculates 
the Kaluza-Klcin (KK) graviton production in supernova cores and the large fraction subsequently 
gravitationally bound around the resulting NS. The predicted decay of the bound KK gravitons to 77 
should contribute to the flux from NSs. Considering 2 to 7 extra dimensions of the same size in the 
context of the HR model, we use Monte Carlo techniques to calculate the expected differential flux of 
gamma-rays arising from these KK gravitons, including the effects of the age of the NS, graviton orbit, 
and absorption of gamma-rays in the magnetosphere of the NS. We compare our Monte Carlo-based 
differential flux to the experimental differential flux using maximum likelihood techniques to obtain 



our limits on LED. Our limits are more restrictive than past EGRET-based optimistic limits that do 
not include these important corrections. Additionally, our limits are more stringent than LHC based 
limits for 3 or fewer LED, and comparable for 4 LED. We conclude that if the effective Planck scale 
is around a TeV, then for 2 or 3 LED the compactification topology must be more complicated than 
a torus. 
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1 Introduction 

In the Standard Model of particle physics, gravity is not unified with the other 3 fundamental forces 
.This is manifested by the hierarchy problem, the fact that the electroweak mass scale Mew ~ 1 TeV 
is many orders of magnitude smaller than the Planck mass scale Mp sa 1.22 x 10 16 TeV [1]. Arkani- 
Hamcd, Dimopoulos, and Dvali (ADD) propose a model of Large Extra Dimensions (LED) as a 
solution for the hierarchy problem. The ADD scenario may be embedded into string theory, which 
allows for the existence of compactified extra dimensions. Due to the presence of n extra dimensions, 
at length scales smaller than the size of the extra dimensions, the gravitational potential between 
test masses has a l/r n+1 dependence; however, on scales larger than the size of the extra dimensions 
the gravitational potential reverts to the ordinary 1/r dependence. For a given n, if all the extra 
dimensions are toroidally compactified, i.e.. have the same size R, the effective Planck mass in the 
(n+4)— dimensional space, Mo, is related to the reduced Planck mass Mp = Mp/s/&K by the relation: 

Ml = R n Mp +2 . (1.1) 

In the ADD model, the hierarchy problem is solved because the presence of LED brings the 
effective Planck mass, M D , to the TeV scale, the truly fundamental scale of gravity. As a consequence, 
the associated Kaluza-Klein gravitons, denoted by Gkk, are massless in the bulk, but they have mass 
on the 3-brane related to their momentum in the bulk (unlike the gravitons of General Relativity). 

According to the ADD model, it is possible to place constraints on extra dimensions by Gkk 
emission from nucleon-nucleon gravibremsstrahlung in type II supernova cores, NN — > NNGkk- 
ADD obtain limits from Supernova (SN) 1987A, assuming pion-exchange mediated by the strong 
force as the dominant process. Hanhart, Reddy and Savage (2001) assume a different process[2]. 
They use nucleon-nucleon gravibremsstrahlung mediated by nucleon-nucleon scattering to obtain the 
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specification 



Fcrmi-LAT 



EGRET 



68% containment PSF (°) at 200 MeV 
Effective Area (cm 2 ) at 200 MeV 
Energy Resolution (%) at 200 MeV 
flux sensitivity (cm" 2 s _1 ) 



2.8 
3000 
13 
6 x 10- 



3.3 
1000 
9.3 
1.3 x 10~ 7 



Table 1: Comparison of performance specifications of Fermi-LAT and EGRET, as relevant for the 
gamma-ray energies considered in this paper[5, 6]. Flux sensitivity is evaluated for a high-latitude 
point source with a E~ 2 spectrum, with 1 year of data, for E > 100 MeV. EGRET effective area is 
quoted for Class A events. 

emission rate for KK gravitons. Furthermore, they indicate that the actual details of the scattering 
process are not important in the soft-radiation limit, where the energy of the outgoing gravitons is 
much less than the energy other incoming nucleons[2]. Then they proceed to obtain limits for n = 2 
and n — 3 extra dimensions from SN1987A. This is based on the argument that the observed neutrino 
luminosity sets an upper bound of 10 19 ergs g^s -1 on the energy loss rate into particles other than 
neutrinos such as Gkk [3]. 

Hannestad and Raffelt (henceforth HR [4]) extend this idea to neutron stars, proposing that if 
the KK gravitons are bound in the gravitational potential of a proto-neutron star as it evolves into 
a neutron star, then the flux of photons from KK graviton decays, Gkk — > 27, could be used to set 
a limit on extra dimensions. They use EGRET results to set limits on LED. However, they do not 
place direct flux limits on the neutron stars not detected by EGRET. Rather, they quote their flux 
upper limit as the 1 yr point-source sensitivity of EGRET for a high latitude point source with a E~ 2 
spectrum (see section 2.1), and derive limits on LED more restrictive than from arguments based 
on KK graviton emissivities from SN 1987A. To obtain upper limits, we follow similar theoretical 
arguments as HR, but we perform a very different analysis, including spectral corrections and upper 
limit spectral analysis with Fermi-LAT data, on 6 gamma-ray faint NS. 

2 Data Analysis 

2.1 Experimental Methods 

The Fermi-LAT is a gamma-ray imaging pair-conversion telescope, consisting of an anti-coincidence 
detector, tracker, calorimeter, and electronics modules. The details of the Fermi-LAT are discussed by 
Atwood et al. [5]. The Fcrmi-LAT features improved performance compared to its predecessor 7-ray 
observatory, EGRET. Some of these specifications, relevant to this study, are compared in Table 1. 

Setting flux limits on sources requires knowledge of background point sources and diffuse emis- 
sion. We make use of the publicly available diffuse models developed by the Fermi-LAT collaboration: 
the Galactic diffuse emission model, glL.iem-v02.fit[7, 8]; and the isotropic model, isotropic-iem_v02.txt [9]. 
The Galactic diffuse model is allowed to vary in a region of interest (ROI) around each source by 
multiplying by a power-law spectral function, as described in [10], effectively making the spectrum 
harder or softer. The scale for the power-law is 100 MeV, and the index is allowed to vary between 
-0.1 and 0.1 (a value of implies no correction to the model). The background point sources are 
modeled by fixing the spectral parameters from the first year Fermi-LAT (1FGL) catalog [10]. 

The data sample consists of a selection of 11 months of all-sky data obtained with the Fermi- 
LAT instrument, using a time interval beginning with the start of survey mode, August 4, 2008, until 
July 4, 2009. This time interval is chosen to be consistent with the 1FGL catalog, so that nearby 
point sources detected with high significance may be modeled appropriately as power- law sources [10] . 
The instrument response function (IRF) chosen is P6- V3-DIFFUSE[ll], as is the case for the 1FGL 
catalog. This IRF specifies a parametrization of effective area, energy resolution, and point spread 
function. We select data from the 1FGL catalog dataset for regions of interest (ROIs) corresponding 
to each source described further in this paper. This dataset excludes events for which the rocking 
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1.24 
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J1136+1551 
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2.13xlO l;l 


J0826+2637 


126.71 


26.62 


196.96 


31.74 


0.53 


0.36 


4.92 


9.64X10 11 



Table 2: Astrophysical properties of neutron star sources analyzed in this work, with sources in 
increasing order of distance. Coordinates, periods, distances, ages, and surface magnetic field 
strengths are obtained from the ATNF Catalog [12]. 

angle is larger than 43°, because of contamination from the Earth's limb due to interactions of cosmic 
rays with Earth's upper atmosphere. For the same reason, for each ROI, events for which the zenith 
angle is larger than 105° are excluded. There is also a good time interval (GTI) cut applied, as 
described in Rcf. [10]. 

2.2 Selection of Neutron Stars 

A query is made on the Australia Telescope National Facility (ATNF) radio pulsar catalog to select NS 
[12]. In order to obtain the best limits, we reject candidate sources that have associations with 7-ray 
sources detected in the 1FGL catalog. To minimize attenuation of the putative KK graviton decay 
signal, as discussed in section IV, we choose NS that satisfy the following criteria: distance d < 0.40 
kpc; surface magnetic field -B SU rf < 5 x 10 13 G; and characteristic age t agc < 2 x 10 8 yr. We take the 
NS ages as the spin-down ages of the pulsars, for consistency over all sources; the corrected ages may 
differ, as discussed in Ref. [13]. However, consideration of the corrected ages hardly affects the limits 
presented here. In addition, the 7-ray sky as viewed by Fermi-LAT is filled with sources near the 
Galactic plane, and diffuse components are also dominant and have large systematic uncertainties at 
low latitudes; we require for Galactic latitude (6), that \b\ > 15° for candidate neutron stars. Applying 
all the above selection criteria to the ATNF catalog, 6 sources remain for analysis, with parameters 
as shown in Table 2. 

2.3 Gamma-ray Spectral Limits 

Fermi-LAT gamma-ray events are selected in a 12° radius ROI centered on each NS source listed in 
Table 2. Given the limitations of the dataset we use and the expected spectral energy distribution 
from gamma rays from trapped KK graviton decay, only gamma-ray events with energies in the range 
100 MeV to 400 MeV are considered. Although desirable, going below 100 McV is not feasible for 
this analysis using P6_V3-DIFFUSE. Before obtaining upper limits for each source, a model for the 
corresponding ROI is developed, inclusive of 1FGL sources and the 2 components of diffuse emission 
(no putative neutron star source is included in this step). 1FGL catalog sources within a 14° radius 
are parametrized as point sources with a power-law spectral energy distribution, with fluxes and 
spectral indices fixed at catalog values; those sources farther than 14° away are not considered. These 
parameters are fixed due to the small gamma-ray energy range (100 MeV-400 McV) that we use, 
which leaves little spectral range to perform accurate fitting. An initial unbinned likelihood fit is 
done in order to determine only the diffuse parameters. For the isotropic diffuse component, the 
parameter to be determined is normalization, while for the Galactic diffuse component, we consider 
the normalization and the spectral index. The analysis, including the diffuse fitting and upper limit 
determination, is performed with the Fermi ScienceTools program pyLikelihood, featuring maximum 
likelihood-based fitting. Version 09-17-00 of the Fermi-LAT ScienceTools is used [14]. For the neutron 
star RX J1856— 3754, a counts map of 100-400 MeV photons in a 10° x 10° region, convolved with a 
Gaussian approximation to the Fermi-LAT PSF, is shown in the left panel of Figure 1. The residual 
counts map (determined from comparison of the counts map to the model-based map), for source 
RX J1856— 3754, is displayed in the right panel of Figure 1. Figure 2 shows a residual counts plot 
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Figure 1: Left: Counts map from data, for source RX J1856— 3754, convolved with a Gaussian 
approximation to the Fermi-LAT PSF, in order to reduce statistical fluctuations without 
dramatically reducing the angular resolution. The colorbar shows counts per pixel. The white 
dashed circle shows the 200 MeV PSF. Right: Residual map, (counts-model) /model, for the same 
source, based on the 1FGL model with the fitted diffuse model. The pixel size is 0.4° for both. 
Green crosses show 1FGL point sources, and the putative 7-ray source is at the center. 
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Figure 2: Residual plot of counts over the 10° square region around source RX J1856— 3754. Black 
horizontal bands indicate the energy range for each point, and black vertical bands represent 
statistical uncertainties, while red vertical bands represent systematic uncertainties of about 10% at 
100 MeV, and decreasing to 5% at 562 MeV[15]. 



versus energy, obtained by integrating over the counts map spatial dependence, and subtracting the 
data counts from model counts and dividing by the model counts. 

A spectral model, which determines the differential flux d^/dE, for each source and number of 



extra dimensions n, is developed in the next section. A significant difference in the data analysis 
technique from Hannestad and Raffelt lies in considering the differential flux rather than the integral 
flux, in determining limits on R; this is a more accurate and optimal method in setting limits in a 
Fermi-LAT analysis, when comparing the data to a pre-defined theoretical distribution. Complete 
details of the theoretical development of the differential flux, as well as analysis methods, can be found 
in [13]. 

3 Calculating the Spectral Model for KK Graviton Decay 7— Rays from 
NS 

In the following subsections, we explain how we calculate the gamma-ray spectrum. Important depar- 
tures from the analysis of HR in forming the theoretical differential flux, d^/dE, to be compared to 
Fermi-LAT observations include: attenuation of the signal due to the age of the neutron star, orbital 
position and velocity of the Gkk, decay of Gkk — > 27, and attenuation of the signal due to mag- 
netic field (which is position and velocity dependent). These features are included via a Monte Carlo 
simulation of about I0 7 Gkk in orbit for each NS source and for n = 2,3, ... ,7 extra dimensions. 

3.1 Theoretical Model 

Following HR, we start with the differential distribution of Gkk created during the proto-neutron 
star core collapse with total energy u> and mass m: 

. ^M NS V NS . (3.1) 

We will make use of notations defined in HR, that we rewrite here for completeness: Q n is the 
total energy loss rate per unit volume into KK gravitons which depends on n; /i — m/oj is the inverse 
of the initial Lorentz factor of the Gkk', At^s — 7.5 s is the time-scale for emission of Gkk during 
the core collapse; and Vns — ^ n R%s 1S ^ ne volume of the proto-neutron star (and current neutron 
star) of radius Rns — 13 km. According to HR, we have: 

^ = Q (RTrn n G n ^)F n (£) , (3.2) 

where R is the extra dimension size, as in eq. (I), T > 30 MeV is the supernova core temperature 
(see Section 5), and f2 n = 27r n / 2 /T(n/2) is the surface of the n-dimensional unit hypersphere, with 
r(„.) as the Gamma function. We also have, 

_ 512 G N an%T^ 

Q ° ~ 5tt3/2 M i/ 2 i 3 - 3 ) 
= 1. 100 x 10 22 McV cm-V 1 (T/30 McV) 7/2 (p/3 x f0 14 g cm" 3 ) 2 (/ X x/0.0075) , (3.4) 

where Newton's constant is Gn = 6.708 x 10~ 33 fic (McV/c 2 )~ 2 , a is the nucleon-nucleon scattering 
cross section of 25 mb, ub is the number density of baryons of 0.16 fm -3 , and M is the isospin- 
averaged nucleon mass of 938 MeV/c 2 . /kk — 0.01 is the estimated fraction of core-collapse energy 
radiated away as Gkk[^G, 17]. In eq. (3.2), the following functions are defined, where q is an integer: 

G q Qi) = ^VT~^ + yM 2 + (3.5) 

f '^ = i'Sn ' (3 ' 6) 

In the previous equation, in writing F(lu/T), we are assuming that the structure function of the 
nuclear medium, in the notation of HR, s(lo/T), of the nuclear medium is unity, which is accurate 



- 5 - 



to first order [4]. An expansion to the next order, (w/T) 2 , would likely shift the expected energy 
distribution of the differential flux to higher energies. Therefore, our assumption of s(u)/T) ~ 1 is in 
the direction of making the associated limits more conservative. Finally, the integral for the case of 
trapped KK gravitons, which make up the initial cloud bound to the NS, is given by: 

N KK , n (t = 0) = 3 / du r 2 dr K /' n . (3.7) 

Jo Jo Ji+u(r) dwd^i 

In eq. (3.7), HR assume that the graviton creation is isotropic at the dimensionless radial 
distance from the neutron star center, r, scaled to the neutron star radius, Rns- The integration over 
r is performed from the proto-neutron star's center to its surface, where r = 1, and the condition 
fi > 1 + U(r) selects the Gkk that are gravitationally trapped. As in HR, we model the neutron 
star's potential as Newtonian: 

GnM ns J (§-±r 2 ) , v<\ 



with U NS = -G N M NS /(R NS c 2 ) = -O.159(Mjvs/l-4M )(13km/ J R iV s). 
The Gkk lifetime is [18] 1 



/100MeV\ 3 
\ m J 



r(m) = 1 x 10 9 yr = K^m" 3 , (3.9) 



where k, = 3.17 x 10 23 MeV 3 s 1 . Assuming an exponential decay of the KK gravitons, the number 
of KK gravitons remaining at time t after the core collapse is given by: 

N K K,n(t) = N K K,n{t = «) exp ( ) (3.10) 

Then, the time derivative (absolute value), of eq. (3.1), is given by: 

A 2 Nkk,u 



d/idw 



= k N KK}n {t) 

Q (RT) n n n At N sV N sKT 2 G n+3 ^)F n+2 (u:/T)e^ (- 1 



r(m) 



3.2 Determining the Differential Flux by Monte Carlo Simulation 

We determine the differential flux according to a Monte Carlo simulation that uses eq. (3.11). We 
calculate the mass distributions and the Lorentz parameters of the decaying gravitons. We then 
determine the momentum and energy distributions of the Gkk, considering the geometry of the 
decays. At the same time, the age of the NS determines the remaining number of gravitons. We also 
consider whether a given gamma ray can escape the NS magnctosphcrc. Finally, this determines the 
differential flux of gamma rays from the NS. 

We carry out the Monte Carlo simulation of the differential flux in the following steps: 

1 ) Sample ui from F„ +2 (cj/T), as in eq. (3.11), for < uj/T < 20 (T = 30 MeV). For oj/T > 20, 
there is negligible contribution from the integral of F n +2(w/T). 

2 ) Sample /i from G„+3(/i), as in eq. (3.11), between fi m in and fi max _. (Note that the sampling steps 
1 & 2 are independent of each other, sec HR.) To simplify the orbit calculation with only a small 
error, we assume that all of the created Gkk start their orbit at the center of the NS (r = 0) 2 . Thus 
we have ^t m i n = 0.807, corresponding to the Gkk escape velocity, and fx max = 0.926, corresponding 



1 This takes into account competing decays to e + e and vv. 

2 We have calculated that this approximation is in the direction of making our limits more conservative. 
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Figure 3: Unit-normalized distributions of KK graviton masses for n = 2, 5, 7, as determined 
according to Monte Carlo simulation. 

to the minimum velocity to reach the neutron star surface, from r — 0. Having determined a value 
of fx, we determine the initial Gkk Lorentz factor 7 = and initial velocity /3 = y/l — p? . Given 
the geometry of the SN explosion, we assume, as do HR, that the Gkk orbits are radial. Using oj 
and [i, from steps 1 and 2, we determine a value for the mass, to = /iuj. Representative distributions 
of to, for different values of n, are shown in Figure 3. 

Since we know the mass at this point, and given the age of the neutron star, t ago , we calculate the 
exponential decay fraction, i'decay = cxp J . Wc sample a real number u uniformly in the 

interval [0,1]: if u > Fdecayj then the event is rejected. 

3 ) Sample the decay vertex ro for a given fi. The probability density function P(r$] fj,), which is shown 
in Figure 4 for two values of [i, is obtained as described in Appendix A. In Figure 5, unit-normalized 
radial profile of decay vertices as a function of radial coordinate, for n = 2, is plotted. 

4 ) Sample the orbit direction isotropically (—1 < cos# < 1, — ir < <j) < ir). This selects an orbit 
direction: 

fo = sin cos 4> x + sin 9 sin cf> y + cos 8 z (3.12) 

where x, y, z are the unit coordinate directions in the NS frame. The z direction is chosen to align 
with the magnetic dipole axis of the NS. At the sampled decay vertex, ro, we then obtain a velocity, 

P' = R NS r/c, (3.13) 
where r = dr/dt\ r=ro is obtained numerically, and the Lorentz factor, 



(3.14) 




(a) fj, = 0.83 (b) fi = 0.90 

Figure 4: Radial probability density functions, P(r;/i), for 2 different values of \i. There are 20 
linearly-spaced bins over the interval [0,r max x P/vs], and the y— axis is the fraction of events per 
bin. 




r x R NS [km] 

Figure 5: Unit-normalized radial distribution of decay vertices as a function of radial coordinate, 
rkm, for n = 2. On the y-axis is plotted (P(r; /z))^, averaged over all /i between ^ m i n and /i max , for 
1 < r < 5.5. 

Using the determined values of m and /?', thus yields p v KK , the 4- momentum in the NS frame 

at the decay vertex. 

5 ) Determine the energy and momentum distribution of one of the two decay photons at decay point 
on the orbit with direction fo- We treat the other decay photon by multiplying the final flux by 
two. Full details of this procedure, as implemented in the Monte Carlo simulation, are given in 
Appendix B. 

6 ) Determine whether the photon pair-produces in the neutron star magnetosphere: this probability 
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is given by P pp (i5 7 , r,py). The probability for photon survival from pair production in the Monte 
Carlo simulation is then taken as P- pp {E 1 ,r,p 1 ) = exp(— t pp ). In Appendix C, we describe the 
computation of t pp . We sample a real number v uniformly on the interval [0,1]: if v > P PP , then 
the event is rejected. 

3.3 The Final Flux Result from the Monte Carlo Simulation 

The distribution defined by the Monte Carlo simulation, dN n /dE 1; is related to the differential flux 
by: 

} =fc„i?"(d/kpc)- 2 ^, (3.15) 



CtEry dE^y 

where the n— dependent constant k n is given as: 

1 ( T \ n 2 

fc » = riT 2 K - -N . n cnrV 1 ™-". (3.16) 

4tt(3.086 x 10 21 ) 3 

and 

N o , n = N KK , n [t = 0)/(RT) n , (3.17) 

where the factor k is related to the decay rate as in eq. (3.9), and the factor T/(hc) is a conversion 
constant, which is numerically 1.52033 x 10 14 m _1 at T = 30 MeV. Values of -ZVo.„ and k n are tabulated 
in Table 3. 

In the computation of dNn/dE^, steps (3) and (7) of Section 3.2 reject events based on the 
decay from the lifetime and the pair production optical depth, respectively. In the case of a zero-age, 
zero-magnetic field neutron star source, the spectrum dN n /dE 1 is normalized to 1. Formally, the 
distribution dN n /dE 7 is defined by: 

dNn^J_dK_ 

dE 1 N ev dE 1 1 ' ' 

where N ev is the number of events in the Monte Carlo simulation. N rem , the number of events 
remaining after the effects of decay and pair production are taken into account, is given by the 
integral: 

/■ 600 MeV ,prr 

^dE^. (3.19) 



dE y 



The upper limit of 600 MeV is determined by the condition lu/T < 20 . The parameter 77, defined as 



400 MeV ^ / ,400 MeV m 



dE~, I ./inn MnV dE. 



100 MeV 



dE-y. (3.20) 

non— attcn 



parameterizes the efficiency with which photons contribute to the spectrum, after signal attenuation 
effects of lifetime and pair production have been taken into account. Values for each source and n arc 
shown in Table 4. 



4 Limits on LED Results 
4.1 Individual Limits 

With all parameters for the ROI fixed, namely catalog sources and diffuse components, as determined 
in Section 2.3, upper limits on R n are determined from spectral fitting based on the method of 
maximum likelihood. Fit values are determined by the MINUIT optimizer [19], and one-sided 95% 
confidence level upper limits are determined by performing a scan of the log-likelihood function in 
each ROI[19]. Statistical parameters of the fit are consistent with non-detection of the KK graviton 

3 The range of gamma-ray energies used to generate dN' n /dE in the Monte Carlo is < E-y < 600 MeV. 
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n 


N .n 


fc„(cm 2 s 1 m n ) 


2 


6.47xl0 4U 


7.126xl0 e 


3 


3.46xl0 41 


5.799xl0 21 


4 


1.94xl0 42 


4.963xl0 ll(5 


5 


7.40 xlO 43 


4.511xl0 51 


6 


7.05xl0 4a 


4.355xlO BB 


7 


4.97xl0 44 


4.452xl0 81 



Table 3: n— dependent constants, as defined in equations (3.17) and (3.16). 



n 


RX J1856-3754 


J0108-1431 


J0953+0755 


J0630-2834 


J1136+1551 


J0826+2637 


2 


0.335 


0.031 


0.221 


0.359 


0.309 


0.332 


3 


0.350 


0.037 


0.249 


0.382 


0.332 


0.360 


4 


0.361 


0.041 


0.276 


0.402 


0.351 


0.385 


5 


0.368 


0.043 


0.302 


0.416 


0.365 


0.406 


6 


0.370 


0.042 


0.325 


0.424 


0.374 


0.419 


7 


0.365 


0.037 


0.334 


0.424 


0.373 


0.423 



Table 4: Table of values of the attenuation parameter 77, defined by eq. (3.20), for the different 

sources analyzed. These attenuation effects can be quite large. HR used only source 

RX J1856-3754 and J0953+0755. These values are calculated for 100 MeV< E 1 < 400 MeV. 




E f [MeV] E v [MeV] 

(a) n = 2 (b) n = 5 

Figure 6: For n = 2 and n = 5, representative distributions of dN n /dEry, according to eq. (3.18), 
for the non-attenuated spectrum and all neutron star sources considered, corrected for magnetic 
pair-production and age attenuation effects. 



decay signal for all NS considered. As a check of this method, we compare upper limits obtained in this 
manner against upper limits computed using profile likelihood implemented by a different method in 
the Fermi-LAT ScienceTools, and find agreement within 10%. Additional systematic checks, due to 
uncertainties in the parameters of the background sources in the ROI, are also performed to verify the 
accuracy of the limits, for source RX J1856— 3754 (the source with the best limits): the agreement in 
the flux upper limits is found to be 15% or better. The flux upper limits for each source and n are 
displayed in Table 5, and the corresponding limits on the LED size R are shown in Table 6. 

4.2 Combined Limits 

We use the following method to combine limits from multiple neutron star sources. A scan over the 
log- likelihood function in each ROI is done with respect to the parameter R n , as shown in Figure 
7. A curve of the change in log- likelihood, |2Alog£|, versus parameter value R n , is generated for 
each source. Then the sum of these curves is taken for all the sources, and the parameter value 
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n 


J1856-3754 


J0108-1431 


J0953+0755 


J1136+1551 


J0630-2834 


J0826+2637 


2 


3.8 


4.3 


5.3 


4.1 


5.8 


6.6 


3 


4.0 


4.4 


5.4 


4.2 


6.8 


8.4 


4 


3.7 


4.2 


6.2 


4.4 


9.4 


9.9 


5 


4.0 


4.1 


6.3 


4.4 


11 


13 


6 


4.1 


4.0 


6.7 


4.2 


14 


15 


7 


4.2 


4.0 


7.8 


3.5 


19 


17 



Table 5: Tabic of 95% C.L. flux upper limits (10 9 cm 2 s 1 ) for the sources analyzed. 



n 


J1856-3754 


J0108-14 


J0953+0755 


J1136+1551 


J0630-2834 


J0826+2637 


2 


9.5 


49 


22 


23 


24 


29 


3 


3.9x10^ 


0.11 


6.7xl0~ 2 


6.9X10" 2 


7.4xl0~ 2 


8.4xl0~ 2 


4 


2.5xl0" 3 


5.4xl0" 3 


3.8xl0- 3 


3.9xl0~ 3 


4.3xl0" 3 


4.8xl0" 3 


5 


5.0xl0~ 4 


9.1xl0~ 4 


7.0xl0~ 4 


7.1xl0~ 4 


8.1xl0~ 4 


8.6xl0~ 4 


6 


1.7xl0~ 4 


2.8xl0~ 4 


2.3xl0~ 4 


2.3xl0~ 4 


2.7xl0~ 4 


2.8xl0~ 4 


7 


8.2xl0~ 5 


1.3xl0~ 4 


1.0xl0~ 4 


1.0xl0~ 4 


1.2xl0~ 4 


1.3xl0~ 4 



Table 6: Table of limits on extra dimensions size R (nm) for the sources analyzed. 




Figure 7: Plot of |2Alog£| versus parameter value of R n for n = 2. The 95% confidence level 
upper limit corresponds to a y— axis value of 2.71, shown by a dashed line. The sum of the curves, 
solid black, is used to obtain the posterior combined limit at the intersection with 2.71. 



corresponding to intersection of that curve with a value of 2.71, corresponding to a one-sided 95% 
confidence level, is quoted as the combined limit value. The results of combining limits on R from 
this method, as well as results from HR, are presented in Table 7. 

4.3 Dependence of LED Limits on Model Parameters 

Dependence of the limits on model parameters, namely T, /kk, and Atjvs, have been evaluated. We 
have determined that the bounds on extra dimensions are quite sensitive to changes in T. Limits 
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R (nm) 


R (nm) 


n 


Combined 


HR 


2 


8.7 


51 


3 


0.037 


0.11 


4 


2.5x10^ 


5.5xl0" 3 


5 


5.0xl0~ 4 


9.1xl0~ 4 


6 


1.7xl0~ 4 


2.8xl0~ 4 


7 


8.2xRT b 


1.2xl0~ 4 



Table 7: 95% CL upper limits on R (nm) for each n, compared to HR2003 EGRET-based results [4]. 



n 


T = 30 MeV 


T = 45 MeV 


2 


9.5 


1.2 


3 


3.9 x 10~ 2 


9.3 x 10" 3 


4 


2.5 x 10- 3 


8.5 x 10~ 4 


5 


5.0 x 10~ 4 


2.1 x 10~ 4 


6 


1.7 x 10~ 4 


8.1 x 10- & 


7 


8.2 x 10~ 5 


4.1 x 10~ 5 



Table 8: A comparison of upper limits on R (nm), evaluated for different values of T, for source 
RX J1856-3754. 



n 




R (nm) 


2 


6.3 


9.5 


3 


8.7 


0.035 


4 


7.4 


2.5 xlO" 3 


5 


5.1 


5.3 xl0~ 4 


6 


9.1 


1.7 xl0~ 4 


7 


9.0 


8.0 xl0~ b 



Table 9: Table of Jkk values and combined upper limits on R (nm) for each value of n, assuming a 
Gaussian prior on Jkk (with mean 0.0075 and sigma 0.00144), as discussed in Sec. 4.4. 



evaluated for source RX J1856— 3754 for T = 30 MeV and a higher value T = 45 MeV, are compared 
in Table 8. The limits on LED are a strong function of temperature. The dependence enters through 
two effects: changing the constant k„ and changing the distribution of gamma-ray energies. The 
limits are affected, since k n ~ y-n-5.5. - m th e r words, by modifying k n , the bounds on LED size 
improve as: 

rp \ — 1— 5.5/n 



In addition, for higher temperatures, the distribution of energies is shifted to higher gamma-ray 
energies. Quantitatively, this increases the integral of the distribution function above 100 MeV, 
fwoMeV^N/dE dE, which tends to improve the bounds. Limits placed on R from source RX 
J1856— 3754 may vary by an order of magnitude, as shown in Table 8. We do not consider lower 
values of T, since according to [3], T = 30 MeV is a conservative lower limit on the SN core temper- 
ature. 

By varying the timescale of core collapse, Atjvs, the limits on R vary as (Aijvs) . Estimates 
for this parameter vary from 5 s to 20 s[20], while we use the value of 7.5 s from HR. Thus, we see 
that the limits depend only weakly on variations of Atjvs and fan- 
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n 


Combined 


CDF 


D0 


LEP 


ATLAS 


CMS 


2 


230 


2.09 


1.40 


1.60 


1.5 


3.2 


3 


16 


1.94 


1.15 


1.20 


1.1 


3.3 


4 


2.5 


1.62 


1.04 


0.94 


1.8 


3.4 


5 


0.67 


1.46 


0.98 


0.77 


2.0 


3.4 


6 


0.25 


1.36 


0.94 


0.66 


2.0 


3.4 


7 


0.11 


1.29 











Table 10: Comparison of 95% CL lower limits on Mo (TeV) with previous astrophysical limits and 
collider limits. Combined limits are obtained in this paper. Collider limits are taken from references 
[22-25]. ATLAS and CMS results arc quoted where A/M D = 1- ATLAS results are quoted with 3.1 
pb _1 of data; CMS results are quoted with 36 pb _1 of data. 

4.4 Effect of Uncertainties on Jxk on the Limits 

Varying the fraction of energy lost into the graviton channel, Jkk, the limits on R vary as f KI ! c . HR 
assumed fan — 0.01, as consistent with diffuse gamma-ray measurements according to EGRET[16, 
17]. However, a more accurate treatment from EGRET low energy diffuse measurements constrains 
Jkk such that 0.005 < Jkk < 0.01. To take this range of values for Jkk into account when computing 
limits, we perform an analysis allowing for a Gaussian prior on the fxK parameter, with a mean of 
0.0075 and a sigma of 0.00144 (as obtained from the variance for a uniform PDF for Jkk between 
0.005 and 0.01). We constrain this parameter to be the same across the 6 ROIs, for each value of 
n. This is possible within the framework of the Fermi-LAT ScienceTools; a similar technique was 
used to constrain dark matter signals from a combined analysis of Milky Way satellites with the 
Fermi-LAT [21]. Limits obtained in this manner are shown in Table 9. 

5 Discussion and Conclusions 

If Md is at a TeV, then for n < 4, the results presented here imply that the compactification topology 
is more complicated than a torus, i.e., all LED having the same size. For flat LED of the same size, 
the lower limits on Mo results are consistent with n > 4. The constraints on LED based on neutron 
star gamma ray emission yield improvements over previously reported neutron star limits, based on 
gamma-ray measurements and combination of individual sources, as shown in Table 10. In addition, 
the results for the n-dimensional Planck mass are much better than collider limits from LEP and 
Tevatron for n < 4, and are comparable or slightly better for n = 4. 

These limits may prove useful, especially for n = 4 case (where the limits are comparable to 
collider results), in the context of constraining phase space in searches for extra dimensions underway 
at the LHC. These results are also more stringent than those reported by short distance gravity 
experiments probing for deviations from the inverse square law. The most sensitive such experiment 
to KK graviton emission presented a result of 37 fim for n = 2 at 95% C.L.[26]; this is several orders 
of magnitude larger than the combined result reported here, of 8.7 nm. 

Casse et al. obtain upper limits significantly better than ours [27] (a factor ~ 20 for n = 2, though 
decreasing approximately as 1/n 2 with increasing n), summing the contribution of all the expected 
NS in the Galactic bulge and comparing to the EGRET data. But it is should be noted that they do 
not account for age nor magnetic field attenuation, while the present analysis shows that both impact 
the photon distribution in a significant way. As these effects are not taken into account in HR, which 
the Casse et al. paper is based on, their upper limits are necessarily underestimated. Furthermore, it 
should be emphasized that the 6 neutron stars analyzed here are chosen at high latitude to avoid the 
large systematic uncertainties involved in modeling the diffuse Galactic background. These are even 
larger in the low energy range that we are interested in here. An analysis of the bulge with the current 
Fermi-LAT instrument response functions and the current Fermi-LAT diffuse models could nominally 
improve the flux upper limits, but at the cost of a much less robust analysis, the systematics of 
which are difficult to evaluate. Within the Fermi-LAT collaboration, a better inner galaxy model is in 
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process, which is necessary before approaching a Galactic Bulge analysis. A Large Extra Dimensions 
analysis of the Galactic Bulge will be the subject of future studies. 



A Appendix: Sampling Decay Vertices from Graviton Trajectories 

In our model, due to the Gkk emission radially outward during the SN core collapse, the Gkk are 
not given any initial angular momentum; thus we assume that the Gkk oscillate on radial paths 
(completely eccentric orbits) through the center of the neutron star. In spherical coordinates (r, 9, <j>), 
this is equivalent to the following constraints: t ) = (j> = 0. The orbital radial distribution, P(r;fj,), is 
defined outside the neutron star by the radial Kepler equation, in which time is given as a function 
of the radial coordinate r[28] 4 : 

t(r) = t k (arcsin (Vkrj - \Jkr{\ - kr) + cij , (A.l) 

where: 

, _ 1 + 1.51^1-7 

fc-Vax- IUns1 , (A.2) 

and: 



_R NS I fc(l-fc) 

tk ~ pck^y i-\u NS \/^- 1 - 6) 



The solution inside the NS (r < 1) is defined as: 



r(t) = ^J=sin(ta), (A.4) 

V\ U NS\ 



where the parameter Q = \/\Uns\c/Rns = 9.13 x 10 3 s _1 . 

Given that t — is the time when the Gkk is created at r = 0, the radial distributions are 
determined by sampling time uniformly between the t = and i max - ^max is given by the time 
to achieve the maximum distance, r max = 1/fc, as in eq. (A. 2). c\ is determined from boundary 
conditions of position and velocity at the surface of the neutron star by solving the full equation of 
motion inside and outside the neutron star. The trajectories for a couple of values of f3 are plotted in 
Figure 8, while the radial distributions for representative values of P{r; fi), are shown in Figure 4. 

B Appendix: Relativistic Decay Kinematics of KK Gravitons 

The energy is given by: 

E 7 = -7'm(l + ^'cos0*), (B.l) 
while the components of the photon momentum, 

P~/=Pxi,~fX + Py'^y + Pz',fz', (B.2) 
in the neutron star frame relative to the direction of the Gkk are given by: 

Px',-y = -msinfl* cos0* (B.3) 
2V, T = — msin#* sin^>* (B-4) 
p^ T = i 7 , m( ( S' + cosr). (B.5) 



4 The radial Kepler equation is not manifestly periodic. However, the full orbit cycle includes an interval over 
< r < r ma x, or a quarter of a cycle. Due to the symmetry of the orbit, our treatment of sampling the decay vertex 
from t = to t(r max ) is sufficient to obtain the full distribution of rg. 
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°o" ' 0.05 oV ' 0.15 ' 0.2 ' ' 0.25 ' 0.3 0.35 ' 0.4 
t[ms] 



Figure 8: Trajectories for /3 = 0.4 and /3 = 0.5. Notice that for larger (3, the KK graviton may 
achieve a farther maximum distance, r max x Rns- This figure shows one quarter of the orbit cycle. 



The z' axis is defined by: z' — fg. 9* is the polar angle between the direction of the Gkk in the lab 
frame (z') and the decay photon in the rest frame of the Gkk, and (jf is the angle in the (x 1 — y') plane. 
In this frame, x' is taken as perpendicular to z' in the z — z' plane and in the direction of increasing 
8, and y' = z' X x' . The coordinate systems used are depicted in Fig. 9. We sample a cosfl* value 
uniformly over the interval [-1,1], and <jf uniformly over the interval [— 7r,7r], given isotropic emission 
of photons in the rest frame of the gravitons. Subsequently, wc obtain momentum components of the 
gamma ray in the frame of the neutron star, p 7 , by rotating back into the frame of the neutron star, 
using the direction of the momentum vector, pkk, as defined by step (5). This is needed for the next 
step. 



C Appendix: Determining Photon Pair Production Optical Depths 

Approximations for the pair-production attenuation are according to the treatment in Refs. [29, 30]. 
The attenuation coefficient depends on the parameter: 



TYl e C B c 

where the critical field is given by 



2 3 

B CI = T^- = 4.414 x 10 13 G, (C.2) 
en 

B± is the magnetic field component of the neutron star perpendicular to the photon's momentum 
vector p 7 , and E 7 is the photon energy. For the magnetic field of the neutron star, we assume a 
static dipole field, 

B(r) = 3 (^-y-" (c.3) 
with dipole moment rh = ^B sur fR^ s z. The attenuation coefficient, a, is given by: 

a{x{E y ,p-y,r)) = -r — -ai(x), (C.4) 

A e £> cl - 

where the reduced attenuation coefficient, ai(x), is expressible as a function of x in terms of a 
modified Bcsscl function of the second kind, with asymptotic limiting expressions for small and large 
values of \ (as plotted in Figure 10): 



X 



Figure 9: The coordinate system, as described in Appendix B. 



i 



« lW = 0.161^3 (A 



) 



0.377e"^,x < 0.1 



(C.5) 



0.597x" 1/3 ,x > 100 



In eq. (C.5), X e = 3.861 x 10~ n cm is the reduced electron Compton wavelength and a.f s is the 
fine structure constant. The asymptotic expressions are used in the Monte Carlo simulation, where 
appropriate, in order to save computer time. 

The optical depth t pp is calculated by path integrating the attenuation coefficient along the 
direction of the photon, from the point of decay, fo =< xq, yo, zq >, out to where r-k m = 7Rns (where 
the field has attenuated to 0.3% of the surface field strength), according to 




(C.6) 



In the preceding equation, s 



given by: 



s 



max — 



xoPi - VoP2 - z a p 3 + J (x pi + yoP2 + z p 3 ) 2 + (7R N s) 2 - xl - yl - z§, 



(C.7) 
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1 



0.1 



0.01 



0.001 



exp(-4/3%) 



10' 



x 10 



10" 



Figure 10: Plot of the reduced attenuation coefficient cti(x) an d limiting asymptotic expressions 
corresponding to cq. (C.5). 



refers to the path length where the photon with direction unit vector * 

P~f = Ttj = Pix + p 2 y + p 3 z (C.8) 
\P-y\ 

is considered to have escaped the magnetosphere. 
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